Transferability of interatomic potentials for silicene

The ability of various interatomic potentials to reproduce the properties of silicene, that is, 2D single-layer silicon, polymorphs was examined. Structural and mechanical properties of flat, low-buckled, trigonal dumbbell, honeycomb dumbbell, and large honeycomb dumbbell silicene phases, were obtained using density functional theory and molecular statics calculations with Tersoff, MEAM, Stillinger–Weber, EDIP, ReaxFF, COMB, and machine-learning-based interatomic potentials. A quantitative systematic comparison and a discussion of the results obtained are reported.


Introduction
We are living in the "silicon age" because of the great importance of elemental silicon to the modern global economy, mainly concerning electronics. Silicon is one of the most extensively investigated materials, and the quality of its production is impressive. However, this applies to bulk silicon. The success of graphene has also sparked interest in other non-carbon 2D materials [1,2]. One of such materials is 2D silicon, called silicene [3,4]. Using first-principles methods with current computer resources enables us to model structures up to about a few hundred atoms. For larger systems, approximate methods are needed, for example, molecular dynamics/statics. For these methods, the quality of the used interatomic potentials (IAPs) is crucial. Because of the importance of silicon, as well as its com-plexity, dozens of potentials have been proposed for it. In the very well-known NIST Interatomic Potentials Repository, there are 27 potentials for silicon (the highest number of potentials provided for a chemical element), the oldest one from 1985 and the latest from 2020 [5,6].
At least five 2D silicon polymorphs have been reported in the literature, that is, flat (FS), low-buckled (LBS) [7], trigonal dumbbell (TDS), honeycomb dumbbell (HDS), and large honeycomb dumbbell (LHDS) silicene [8]. There are still doubts about their dynamic stability. For example, for a flat phase, the negative ZO phonon mode could be removed by the selection of an appropriate substrate [3,9]. The ability of poten-tials for 3D silicon to reproduce 2D silicon is poorly studied. There are several papers in which the quality of potentials for 3D silicon has been assessed [10][11][12], but not for silicene. The intention of this work is first to determine the structural and mechanical properties of 2D silicon using the first-principles method and then to test the ability of different interatomic potentials to reproduce these properties.

Methods
Analyzing the available literature concerning all phases of single-layer Si, it is not feasible to find all structural, mechanical, and phonon data obtained in one consistent way. The availability of experimental data is actually limited to the silicene grown on supports, a pristine free-standing single-layer sheet of silicene has not yet been discovered [4,13]. Therefore, we must use ab initio calculations. Unfortunately, also ab initio calculations, most often DFT, differ in the calculation methodology, that is, they use different functional bases and different pseudopotentials or exchange-correlation functionals. Also, parameters such as cohesive energy and elastic constants are poorly accessible. For this reason, structural and mechanical data, that is, lattice parameters, average cohesive energy, average bond length, average height, 2D elastic constants, as well as phonon data are determined here using a single consistent firstprinciples approach as described in the next section "Ab initio calculations". These data were further considered as reference data and marked as "value DFT ". Then, the same data were determined as described in section "Molecular calculations" using the analyzed molecular potentials from subsection "Interatomic potentials" and are marked as "value potential ". Having both data, we can simply define the mean absolute percentage error (MAPE): (1) which enables us to quantify the potentials under examination.
For 2D materials, directional 2D Young's moduli, 2D Poisson's ratios, (3) and the 2D shear modulus, (4) are often used instead of elastic constants C ij . Because of the symmetry of hexagonal lattices, these reduce to one 2D Young's modulus E and one 2D Poisson's ratio ν [14].

Ab initio calculations
The ab initio calculation methodology here is closely analogous to that used in [15]. Hence, its description is also very similar, that is, density functional theory (DFT) [16,17], ABINIT plane-wave approximation code [18,19], local density approximation (LDA) [20,21] as an exchange-correlation functional, and optimized norm-conserving Vanderbilt pseudopotential [22] (ONCVPP) are similar. Cut-off energy and electron configuration of Si were used in the DFT calculations according to the pseudopotential and Gaussian smearing scheme with tsmear (Ha) = 0.02. To generate k-points grids, kptrlen was set to 43.0. Since the 2D structures were analyzed in the z direction, a vacuum of 20 Å was applied. The initial data for the five structures analyzed were deduced from [7] and [8]. The structures were then carefully relaxed with full optimization of cell geometry and atomic coordinates [15].
The average cohesive energy E c (eV/atom) was computed as the difference in the total energy of a given relaxed earlier structure and that of its individual atoms placed in a cubic box of sufficient size. The theoretical ground state, T = 0 K, and the elastic constants, C ij , of all previously optimized structures were computed using the metric tensor formulation of strain in the density functional perturbation theory (DFPT) [23]. The mechanical stability of the analyzed structures was verified by calculating the so-called Kelvin moduli [24,25]. To calculate the phonons, the DFPT implemented in ABINIT [18,19] [26] of the analyzed structures were then used to identify their dynamical stability [27], complementary to the mechanical stability.
As for DFT calculations, the structures here were fully prerelaxed with the conjugate gradient (CG) algorithm and, then, the elastic constants, C ij , were calculated for them using the stress-strain method with the maximum strain magnitude set to 10 −6 [30,31]. In the z direction, a vacuum was set to 20 Å.
To measure the performance of the analyzed interatomic potentials, a series of molecular dynamics (MD) simulations (200 atoms and 10000 timesteps, NVE ensemble) and LAMMPS's built-in function timesteps/s were used. The results were then normalized relative to the longest run time.

Interatomic potentials
The parameterizations of the potentials listed below were obtained from the NIST Interatomic Potentials Repository and/or from LAMMPS code sources.
1. Tersoff1988 [33]: the original Tersoff potential for silicon (it is important to remember that this paper proposed a form of potential rather than a specific parametrization for silicon) 2. Tersoff2007 [34]: the Tersoff potential fitted to the elastic constants of diamond silicon 3. Tersoff2017 [35]: newer, better optimized the Tersoff potential for silicon 4. MEAM2007 [36]: the semi-empirical interatomic potential for silicon based on the modified embedded atom method (MEAM) formalism 5. MEAM2011 [37]: the spline-based modified embeddedatom method (MEAM) potential for Si fitted to silicon interstitials 6. SW1985 [38]: the Stillinger-Weber (SW) potential fitted to solid and liquid forms of Si 7. SW2014 [39]: the Stillinger-Weber (SW) potential fitted to phonon dispersion curves of a single-layer Si sheet 8. EDIP [40]: the environment-dependent interatomic potential (EDIP) fitted to various bulk phases and defect structures of Si 9. ReaxFF [41]: the reactive force-field (ReaxFF) fitted to a training set of DFT data that pertain to Si/Ge/H bonding environments 10. COMB [42]: the charge optimized many-body (COMB) potential fitted to a pure silicon and five polymorphs of silicon dioxide 11. SNAP [43]: the machine-learning-based (ML-IAP) linear variant of spectral neighbor analysis potential (SNAP) fitted to total energies and interatomic forces in groundstate Si, strained structures, and slab structures obtained from DFT calculations 12. qSNAP [43]: the machine-learning-based (ML-IAP) quadratic variant of spectral neighbor analysis potential (qSNAP) fitted to total energies and interatomic forces in ground-state Si, strained structures, and slab structures obtained from DFT calculations 13. SO(3) [44]: the machine-learning-based (ML-IAP) variant of the SO(3) smooth power spectrum potential (SO (3)) fitted to the ground-state of the crystalline silicon structure, strained structures, slabs, vacancy, and liquid configurations from DFT simulations 14. ACE [45]: the machine-learning-based (ML-IAP) variant of the atomic cluster expansion potential (ACE) fitted to a wide range of properties of 3D silicon determined from DFT calculations

Results and Discussion
Structural and mechanical properties  Figure 1. Additionally, the crystallographic data for them are stored in crystallographic information files (CIFs) in Supporting Information Files 1-5.
The results of first-principles calculations show that all silicene phases have hexagonal symmetry. The symmetry characteristics of a structure determine the symmetry of its physical properties (cf. Neumann's Principle and Curie laws) [30,46]. For 2D linear hyperelastic materials, there are four classes of symmetry [25], and hexagonal symmetry implies isotropy of the stiffness tensor, that is, there are only two distinct elastic constants and they satisfy, in Voigt notation, such conditions that C 11 = C 22 , C 33 = (C 11 − C 12 )/2.
Structural and mechanical characteristics that were determined from DFT computations, namely lattice parameters, average cohesive energy, average bond length, average height, 2D elastic constants, 2D Young's modulus, Poisson's ratio, and 2D Kelvin moduli, of the five silicene polymorphs analyzed are gathered in Table 1. Since we are analyzing free-standing silicene here, which has not yet been observed in experiments, we compare the results of the calculations with those of other authors. We find that the lattice constants, average bond length, average height and cohesive energy agree at the DFT level of accuracy with other calculations. Mechanical properties of silicene are available in the literature only for the LBS phase and are limited to 2D Young's modulus and Poisson's ratio only. These quantities are also in reasonable agreement with the present results. It is worth noting that all calculated 2D Kelvin moduli for all silicene phases are positive, which results in mechanical stability [25].
The phonon spectra along the Γ-M-K-Γ path for the five silicene polymorphs is depicted in Figure 2. The analysis of the computed curves shows that the phases TDS, LBS, HDS, and LHDS are not only mechanically but also dynamically stable, that is, all phonon modes have positive frequencies anywhere. The FS phase is mechanically stable, but can be dynamically unstable, that is, the optical ZO phonon mode has a negative frequency. Other authors have also observed similar FS phase behavior [4,13]. However, since silicene is not a free-standing structure in nature, the selection of a proper substrate may dampen the out-of-plane vibration mode, and flat silicene may be produced [3]. Hence, it was decided to include also this phase in the molecular calculations.

Performance of interatomic potentials
Computations were carried out using molecular statics and the fourteen interatomic potentials for silicon, (Tersoff (×3), MEAM (×2), Stillinger-Weber (×2), EDIP, ReaxFF, COMB and machine-learning-based (ML-IAP (×4)), enumerated in section "Interatomic potentials". Twelve structural and mechanical properties, that is, lattice parameters a and b, average cohesive energy E c , average bond length d, average height h, 2D elastic constants C ij , and 2D Kelvin moduli K i , are collected in Table 2 for flat silicene (FS), in Table 3 for low-buckled silicene (LBS), in Table 4 for trigonal dumbbell silicene (TDS), in Table 5 for honeycomb dumbbell silicene (HDS), and in Table 6 for large honeycomb dumbbell silicene (LHDS), respectively. The aforementioned results, for each of the five silicene phases, were then compared with those from DFT calculations using the mean absolute percentage error (MAPE) defined in Equation 1. Let us briefly analyze the results for each phase. For the FS phase, the most accurate is the MEAM2011 potential as it has the lowest MAPE , see Table 2. For the LBS phase, it is the Tersoff2107 potential, see Table 3. For the TDS phase, it is the ReaxFF potential, see Table 4. For the HDS phase, it is the ReaxFF potential, see Table 5, and finally, for the LHDS phase, it is again the ReaxFF potential, see Table 6. Now let us take a summary look. Seven of the analyzed fourteen potentials, namely Tersoff2007, Tersoff2017, SW1985, SW2014, ReaxFF, SNAP, and ACE, are able to correctly reproduce the structural properties of the five polymorphs of silicene, see Table 3, Table 5, and Table 6. Two potentials, ReaxFF and MEAM2011, give the best quantitative performance measured by the total mean absolute percentage error (MAPE), see Table 6. Regarding the cost of calculations in terms of relative performance measured as normalized timesteps per second in         molecular dynamics (MD) simulations, the EDIP and Stillinger-Weber potentials are the fastest, about five times faster than the MEAM and Tersoff potentials, about 100 times faster than ReaxFF and COMB, and up to 2000 times faster than the ML-IAP(ACE) potential, see Table 6. It is also worth noting, that the machine-learning-based (ML-IAP) interatomic potentials, according to the methodology used, are not superior to classical potentials in terms of performance (MAPE) and are instead even three orders of magnitude more computationally expensive, see Table 6.

Conclusion
A systematic quantitative comparative study of various silicon interatomic potentials for reproducing the properties of five silicene (2D silicon) polymorphs was shown. In order to compare the fourteen potentials listed in section "Interatomic potentials", the structural and mechanical properties of flat (FS), low-buckled (LBS), trigonal dumbbell (TDS), honeycomb dumbbell (HDS), and large honeycomb dumbbell (LHDS) silicene ( Figure 1) obtained from density functional theory (DFT) and molecular statics (MS) computations were used. Computational cost and performance of the analyzed potentials were compared.
Considering the performance and the cost of calculations, the classical potentials of Tersoff, SW, and MEAM types seem to be the best choice here. Although data for silicene polytypes were not used in the optimization of these potentials, they were able to reproduce their properties well. This is a consequence of the fact that they are based on physics, have a natural extrapolation ability, and do not just interpolate data.
I hope that the findings presented here will help other researchers in selecting the suitable potentials for their purposes and will be a hint to parameterize new potentials for silicene.

Supporting Information
Crystallographic Information Files (CIF) for polymorphs of silicene (created by qAgate (Opensource software to post-process ABINIT) and then converted to P1 setting by Bilbao Crystallographic Server [48,49]